import holoviews as hv
hv.extension('bokeh')
hv.opts.defaults(hv.opts.Curve(width=500), 
                 hv.opts.Histogram(width=500))

import numpy as np
import scipy.stats

#Sin curve example
import itertools
import matplotlib.pyplot as plt
from sklearn.mixture import GaussianMixture
from scipy import linalg
import matplotlib as mpl

print(f"NumPy versión: {np.__version__}")
print(f"Scipy versión: {scipy.__version__}")
NumPy versión: 1.21.5
Scipy versión: 1.9.1

20. Modelos de Mezcla de Gaussianas#

20.1. Modelos de variable latente (LVM)#

Métodos no supervisados: Contrario a los modelos supervisados que vimos antes (regresión lineal, regresión logística), en el escenario no supervisado, no se tienen etiquetas o variable objetivo a ajustar.

El objetivo general de los métodos no-supervisados es revelar estructura en los datos y así facilitar su entendimiento.

Un ejemplo de lo anterior es la identificación de una causa oculta o variable latente. Este es el dominio de los modelos de variable latente (Latent Variable Models, LVM).

En un LVM se tiene la siguiente estructura:

  • Variable observada (efecto): Corresponde a nuestros datos.

  • Variable latente (causa): Una variable que busca explicar la variable observada.

Y en general el objetivo es: inferir la variable latente dada la variable observada. Para lo anterior asumimos alguna propiedad para la variable latente, por ejemplo:

  • La dimensión de la variable latente es menor que la de la variable observada: Métodos de reducción de dimensionalidad.

  • La variable latente es continua (e.g., PCA, FA).

  • Las variables latentes son independientes (la variable observada si puede estar correlacionada).

En esta lección nos enfocaremos en un tipo de LVM con variable latente discreta, llamado modelo de mezcla de gaussiana, con dos aplicaciones:

La estimación de densidad y generación de datos

Clustering (agrupamiento): Agrupar datos que sean “similares” en conjuntos (clusters)

20.2. Modelo de mezcla de gaussiana (GMM)#

20.2.1. Modelo de mezcla#

Un modelo de mezcla (mixture model) general es un LVM con variable latente discreta donde cada dato (variable observada) se representa como:

\[ x_i \in \mathbb{R}^D \text{ with } i=1,\ldots, N \]

y tiene asociada una categoría discreta (variable latente):

\[ z_i \in {1, 2, \ldots, K} \]

Podemos modelar el prior de \(z_i\) como una distribución categórica:

\[ f(z_i = k) = \pi_k \]

donde \(\pi_k \in [0, 1]\) son los coeficientes de mezcla y \(\sum_{k=1}^K \pi_k=1\).

Con lo anterior podemos escribir la verosimilitud del modelo de mezcla para cada dato como:

(20.1)#\[\begin{split} \begin{align} L(\theta) = f(x_i ; \theta) &= \sum_{k=1}^K f(z_i=k; \theta_k) f(x_i|z_i = k; \theta_k) \\ &= \sum_{k=1}^K \pi_k f(x_i|z_i = k; \theta_k) \end{align} \end{split}\]

Importante

Estamos asumiendo que existen \(K\) componentes (conjuntos, clusters).

Para continuar debemos especificar una familia para \(f(x_i|z_i = k; \theta_k)\).

Advertencia

En esta lección, las notaciones se simplifican sin diferenciar escaladores, vectores y matrices con negrita, con el fin de centrarse en las ideas.

20.2.2. Modelo de mezcla de gaussiana (GMM)#

Para el modelo de mezcla de gaussiana (Gaussian Mixture Model, GMM), asumiremos que cada componente de la mezcla sigue una distribución Gaussiana/normal.

Con lo anterior la Ecuación (20.1) para cada observación queda como:

\[ L(\theta) = f(x_i ; \theta) = \sum_{k=1}^K \pi_k \mathcal{N}(x_i; \mu_k, \Sigma_k) \]

donde cada componente \(k\) tiene tres parámetros:

  • Media: \(\mu_k\). Representa la ubicación del centro del componente.

  • Covarianza: \(\Sigma_k\): Representa el tamaño y la orientación del componente.

  • Concentración: \(\pi_k\): Representa la probabilidad del componente.

Luego ajustar el GMM se reduce a encontrar los mejores valores para estos parámetros para todos los componentes dado los datos

Nota

El GMM es un modeloge generativo, dado un valor para los parámetros podemos muestrar del GMM, generando datos nuevos.

Exploremos la influencia de estos parámetros a continuación con un ejemplo unidimensional y luego otro bidimensional. Asocie los valores de los parámetros con las gaussianas que observa en las figuras.

from utils import create_mixture

data, labels = create_mixture(pis=[0.4, 0.5, 0.1], 
                              mus=[[-1], [2], [2.5]], 
                              Sigmas=[0.1, 0.5, 0.1])

edges, bins = np.histogram(data, bins=50, density=True)
Hide code cell source
hv.Histogram((edges, bins), kdims='data', vdims='Density')
data, labels = create_mixture(pis=[0.9, 0.1], 
                              mus=[[-1, -1], [1, 1]], 
                              Sigmas=[0.5*np.eye(2), [[1, -0.9],[-0.9, 1]]])

bins, edgesx, edgesy = np.histogram2d(data[:, 0], data[:, 1], bins=30, density=True)
Hide code cell source
hv.Scatter((data[:, 0], data[:, 1])) +  hv.Image((edgesx, edgesy, bins))

20.2.3. Probabilidad a posteriori del GMM#

Asumiendo que conocemos los parámetros del GMM, una aplicación relevante es la de inferir el componente más probable para un dato en particular,

Esto está dado por la probabilidad de \(z\) dado \(x\).

Utilizando la regla de Bayes y los supuestos previos podemos escribir:

(20.2)#\[\begin{split} \begin{split} r_{ik} = f(z_i = k|x_i; \theta) &= \frac{ f(x_i|z_i = k; \theta) f(z_i=k; \theta)}{f(x_i)} \\ &= \frac{ f(x_i|z_i = k; \theta) f(z_i=k; \theta)}{ \sum_{k=1}^K f(x_i|z_i = k; \theta) f(z_i=k; \theta) } \\ &= \frac{ \mathcal{N}(x_i; \mu_k, \Sigma_k) \pi_k}{ \sum_{k=1}^K \mathcal{N}(x_i; \mu_k, \Sigma_k) \pi_k }, \end{split} \end{split}\]

este posterior también suele llamarse responsabilidad \(r_{ik}\) del componente \(k\) para el dato \(i\). Esto es lo que se conoce como una asignación suave (probabilidad de pertenecer a cada componente)

En algunos casos sólo necesitamos una asignación fuerte, es decir descubrir cual es el componente más probable:

\[ c_{i} = \text{arg}\max_k r_{ik} = \text{arg}\max_k \left(\log \mathcal{N}(x_i;\mu_k, \Sigma_k) + \log \pi_k\right) \]

En cuyo caso podemos omitir la evidencia, es decir el denominador de la Ecuación (20.2).

La siguiente figura muestra la diferencia graficamente la diferencia entre una asignación suave y dura, en este caso para dos componentes.

def jointGMM(data, pi, mu, Sigma):
    # Assume 2 dimensional data
    data_centered = data - mu
    norm = pi/(2*np.pi*np.sqrt(np.linalg.det(Sigma)))
    xSx = np.multiply(data_centered, np.linalg.solve(Sigma, data_centered.T).T)
    return norm*np.exp(-0.5*np.sum(xSx, axis=1))

pis = [0.6, 0.4]; mus=[[-1, -1], [1, 1]]; 
Sigmas=[0.5*np.eye(2), [[1, -0.3],[-0.3, 1]]]
data, labels = create_mixture(pis, mus, Sigmas, rseed=0)

posterior = np.zeros(shape=(len(data), len(pis)))
for k in range(len(pis)):
    posterior[:, k] = jointGMM(data, pis[k], mus[k], Sigmas[k])
posterior = posterior/np.sum(posterior, axis=1)[:, np.newaxis]
Hide code cell source
p1 = hv.Scatter((data[:, 0], data[:, 1], posterior[:, 1]), 
                vdims=['y', 'z'], label='Soft assignment').opts(color='z')
p2 = hv.Scatter((data[:, 0], data[:, 1], np.argmax(posterior, axis=1)), 
                vdims=['y', 'z'], label='Hard assignment').opts(color='z')
hv.Layout([p1, p2]).cols(2).opts(hv.opts.Scatter(cmap='RdBu'))

20.2.4. Ajustando el GMM#

Ajustar se refiere a encontrar los estimadores puntuales óptimos para el modelo estadístico.

Anteriormente aprendimos sobre el procedimiento de MLE. El logaritmo de la verosimilitud en este caso es:

\[ \log \mathcal L(\pi,\mu,\Sigma) = \sum_{i=1}^N \log \sum_{k=1}^K \pi_k \mathcal{N}(x_i;\mu_k, \Sigma_k) \]

Advertencia

Con la excepción de \(K=1\), no se puede encontrar una expresión analítica para los valores óptimos de los parámetros.

En general cuando no tenemos datos completos, es decir cuando tenemos variables ocultas (latent) o perdidas (missing), el posterior no es sencillo de factorizar.

En este caso particular tenemos una variable latente \(z_i\), la cual necesitamos marginalizar (sumar sobre \(k\)) para calcular la verosimilitud. En lo que sigue veremos un método de estimación aproximada conocido como Expectation Maximization (EM).

Algunos otros problemas del MLE/MAP para el caso de GMM:

  • No-convexo: muchos óptimos locales (Murphy 11.3.2)

  • Muy difícil de calcular, es NP-hard

  • Unidentifiability: Puede existir \(K!\) configuraciones con exactamente la misma solución (label switching).

20.3. Método de Expectation Maximization (EM)#

La solución de MLE requiere que podemos observar completamente los datos. Esto significa que

Conocemos \(z_i\) para cada \(x_i\), es decir que conocemos a componente pertenece \(x_i\)

Antes en el ejemplo de GMM, tenemos log-verosimilitud parcialmente observada (es decir, sin saber el valor de la variable latente \(z_i\)):

\[ \log \mathcal L(\theta) = \log \prod_{i=1}^N f(x_i; \theta) = \sum_{i=1}^N \log f(x_i; \theta) = \sum_{i=1}^N \log \sum_{k=1}^K f(z_i=k;\theta_k) f(x_i|z_i=k; \theta_k) \]

La log-verosimilitud completamente observada (fully-observed, FO) es:

\[ \log \mathcal L_{\text{FO}}(\theta) = \log \prod_{i=1}^N f(x_i, z_i; \theta) = \sum_{i=1}^N \log f(x_i, z_i; \theta) = \sum_{i=1}^N \log f(z_i;\theta) f(x_i|z_i; \theta) \]

sin embargo en la práctica no conocemos \(z\), por lo tanto la marginalizamos calculando el valor esperado dado el posterior de \(z\):

\[\begin{split} \begin{split} Q(\theta | \theta^{\text{old}}) &= \mathbb{E}_{f(z|x; \theta^{\text{old}})} \left[\log \mathcal L_{\text{FO}}(\theta) \right] \\ &= \sum_{k=1}^K \sum_{i=1}^N f(k|x_i; \theta^{\text{old}}_k) \log f(k;\theta_k) f(x_i|k; \theta_k) \end{split} \end{split}\]

que se conoce como la función auxiliar. Luego de esto actualizamos los parámetros del modelo maximizando:

\[ \theta^{\text{new}} = \text{arg}\max_\theta Q(\theta | \theta^{\text{old}}), \]

y finalmente asignamos

\[ \theta^{\text{old}} \leftarrow \theta^{\text{new}} \]

repitiendo este procedimiento hasta tener convergencia.

Este procedimiento iterativo:

  • Paso E: Estimar el valor esperado de la verosimilitud dado los parámetros actuales

  • Paso M: Maximizar el estimador y actualizar los parámetros

se llama Expectation Maximization (EM). EM es un algoritmo general:

  • puede utilizarse para resolver otros problemas sin solución analítica o con variables no observadas.

  • tiene muchas variantes: Incremental, Variacional, Monte-Carlo (Murphy 11.4.8 and 11.4.9)

20.3.1. Actualizaciones de EM para el caso de GMM#

La función auxiliar para GMM sería:

\[\begin{split} \begin{split} Q(\theta | \theta^{\text{old}}) &= \sum_{k=1}^K \sum_{i=1}^N f(k|x_i; \theta^{\text{old}}_k) \log f(k; \theta_k) f(x_i|k; \theta_k) \\ &= \sum_{k=1}^K \sum_{i=1}^N r_{ik}^{\text{old}} \left(\log \pi_k + \log \mathcal{N}(x_i; \mu_k, \Sigma_k) \right) \end{split} \end{split}\]

donde \(r_{ik}^{\text{old}}\) es la responsabilidad dado los parámetros de la iteración previa. Luego maximizamos \(Q\) para actualizar los parámetros.

El estimador para las medias está dado por:

\[ \mu_{k}^{\text{new}} = \frac{ \sum_{i=1}^N r_{ik}^{\text{old}} x_i}{\sum_{i=1}^N r_{ik}^{\text{old}}} \]

El cual se obtiene de derivar \(Q\) e igualarlo a cero.

Demostración

La función de densidad de una distribución normal multivariante con D dimensión es:

\[ f(x; \mu, \Sigma) = \frac{1}{(2\pi)^{\frac{D}{2}}|\Sigma|^{\frac{1}{2}}} exp\left\{-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu) \right\} \]

Hacemos la derivación:

\[\begin{split} \begin{align} \frac{d Q}{d\mu_{k'}} &= \frac{d}{d\mu_{k'}} \sum_{k=1}^K \sum_{i=1}^N r_{ik}^{\text{old}} \left(\log \pi_k + \log \mathcal{N}(x_i;\mu_k, \Sigma_k) \right) \nonumber \\ &= \sum_{i=1}^N r_{ik'}^{\text{old}} \frac{d}{d\mu_{k'}} \log \mathcal{N}(x_i;\mu_{k'}, \Sigma_{k'})\nonumber \\ &= \sum_{i=1}^N r_{ik'}^{\text{old}} \frac{d}{d\mu_{k'}} \left( -\frac{1}{2} (x_i - \mu_{k'})^T \Sigma_{k'}^{-1} (x_i - \mu_{k'}) -\frac{D}{2} \log(2\pi) -\frac{1}{2} \log(|\Sigma_{k'}|) \right)\nonumber \\ &= -\frac{1}{2} \sum_{i=1}^N r_{ik'}^{\text{old}} \frac{d}{d\mu_{k'}} (x_i - \mu_{k'})^T \Sigma_{k'}^{-1} (x_i - \mu_{k'}) \nonumber \\ &= \sum_{i=1}^N r_{ik'}^{\text{old}} \Sigma_{k'}^{-1} (x_i - \mu_{k'}) = 0 \nonumber \end{align} \end{split}\]

Siguiendo el mismo procedimiento para la covarianza se llega a:

\[\begin{split} \begin{align} \Sigma_k^{\text{new}} &= \frac{1}{\sum_{i=1}^N r_{ik}^{\text{old}}} \sum_{i=1}^N r_{ik}^{\text{old}} (x_i - \mu_k) (x_i - \mu_k)^T \nonumber \\ &= -\mu_k \mu_k^T + \frac{\sum_{i=1}^N r_{ik}^{\text{old}} x_i x_i^T}{\sum_{i=1}^N r_{ik}^{\text{old}}} \nonumber \end{align} \end{split}\]

Finalmente, incorporando la restricción \(\sum_{k=1}^K \pi_k = 1\) se llega a:

\[ \pi_k^{\text{new}} = \frac{\sum_{i=1}^N r_{ik}^{\text{old}}}{N} \]

20.3.2. Resumen del algoritmo de EM para ajustar GMM#

  • (1) Inicializar (aleatoriamente) los centros de los componentes.

  • (2) Calcular la responsabilidad de los componentes para cada dato con:

\[ r_{ik}^{\text{old}} = \frac{ \mathcal{N}(x_i; \mu_k^{\text{old}}, \Sigma_k^{\text{old}}) \pi_k^{\text{old}}}{ \sum_{k=1}^K \mathcal{N}(x_i; \mu_k^{\text{old}}, \Sigma_k^{\text{old}}) \pi_k^{\text{old}} } \]
  • (3) Calcular los nuevos valores de los parámetros con:

\[ \mu_{k}^{\text{new}} = \frac{\sum_{i=1}^N r_{ik}^{\text{old}} x_i}{\sum_{i=1}^N r_{ik}^{\text{old}}} \]
\[ \Sigma_k^{\text{new}} = -\mu_k \mu_k^T + \frac{\sum_{i=1}^N r_{ik}^{\text{old}} x_i x_i^T}{\sum_{i=1}^N r_{ik}^{\text{old}}} \]

y

\[ \pi_k^{\text{new}} = \frac{\sum_{i=1}^N r_{ik}^{\text{old}}}{N} \]
  • (4) Actualizar los parámetros con \(\theta^{\text{old}} \leftarrow \theta^{\text{new}}\).

  • (5) Si se cumple el criterio de convergencia entonces terminar, de lo contrario volver a 2.

Aquí es una visualización del proceso con iteraciones (L) paso a paso (ref: Bishop Fig. 9.8):

../../_images/gmm_bishop.png

20.3.3. Consideraciones prácticas#

Hay varias consideraciones prácticas para EM (inicialización, Convergencia) que se aplican al GMM. También hay una consideración de número de componentes de GMM en la prática.

Inicialización

El procedimiento de EM garantiza convergencia a un punto estacionario. Pero está sujecto a un óptimo local. Se puede obtener una mejor solución con las siguiente heurísticas:

  • Realizar varias inicializaciones (e.g., aleatorias) de los párametros de diferentes áreas y elegir el modelo con el mejor rendimiento.

  • Inicializar con la solución de un modelo más simple (por ejemplo k-means para GMM). Esto puede evitar errores catastróficos pero también puede limitar la exploración. (Ojo: k-means también se puede mejorar mediante varias inicializaciones.)

Convergencia

EM es un procedimiento iterativo. ¿Cuántas veces repetimos el procedimiento? Una heurística típica es detener el algoritmo si el logaritmo de la verosimilitud crece menos que un cierto valor porcentual con respecto al paso anterior (estancamiento).

Número de componentes de GMM

El número de componentes \(K\) es un hiperparámetro que se escoje a priori por el usuario. ¿Cómo saber que valor de \(K\) es adecuado?

Una opción es utilizar métricas que consideran tanto error como complejidad como por ejemplo

  • Criterio de información de Akaike (Akaike Information Criterion, AIC)

  • Criterio de información bayesiana (Bayesian Information Criterion, BIC)

\[ \text{BIC} = K \log N - 2 \log \mathcal{L} \]

Una tercera opción es utilizar el coeficiente de silueata.

20.4. GMM con la librería scikit-learn#

Del módulo sklearn.mixture tenemos la siguiente clase que implementa una mezcla de Gaussianas:

GaussianMixture(n_components=1, # Number of components
                covariance_type='full', # Shape of the covariance
                tol=0.001, # termination tolerance of EM
                max_iter=100, # Max number of EM iterations
                n_init=1, # Number of random initializations
                init_params='kmeans', # Initialization for the responsabilities (kmeans or random)
                ...
               )

donde el argumento covariance_type puede ser full, tied, diag o spherical.

  • full: Las covarianzas completas significa que no hay restricciones en las matrices de covarianza.

  • tied: Las covarianzas enlazadas significa que todos los componentes tienen la misma matriz de covarianza (comparten parámetros).

  • diag: En una matriz de varianza diagonal, las dimensiones no tienen correlación (es decir, los componentes no pueden tener rotaciones), pero la varianza de cada dimensión puede ser diferente.

  • spherical: Las dimensiones no tienen correlación (es decir, los componentes no pueden tener rotaciones), y todas las dimensiones tienen igual varianza. Es decir, para cada componente \(k\), \(\Sigma_k = \sigma_k^2 \text{I}\), \(\sigma_k \in \mathbb{R}_{+}\).

../../_images/tipo_var_covar.png

Probemos el módulo utilizando datos sintéticos:

data, labels = create_mixture(pis=[0.33, 0.33, 0.34], 
                              mus=[[-3, 0], [0, 0], [3 ,0]], 
                              Sigmas=[[[1, 0.9],[0.9, 1]], 
                                      [[1, -0.9],[-0.9, 1]], 
                                      [[1, 0.9],[0.9, 1]]])
Hide code cell source
hv.Scatter((data[:, 0], data[:, 1])).opts(width=500)

20.5. Uso 1: Clustering#

Usamos los datos sintéticos arriba suponiendo que no sabemos cómo se generan, hacemos una tarea de clustering ahora!

Utilizaremos el BIC para encontrar el “mejor” número de clusters.

  • Para ajustar el modelo utilizamos el método fit

  • Para calcular el BIC utilizamos el método bic

from sklearn.mixture import GaussianMixture

maxK = 20; 
cov = 'full' 
metrics = np.zeros(shape=(maxK, ))
for i, K in enumerate(range(1, maxK+1)):
    gmm = GaussianMixture(n_components=K, covariance_type=cov, n_init=3)
    gmm.fit(data)
    metrics[i] = gmm.bic(data)
Hide code cell source
hv.Scatter((range(1, maxK+1), metrics), 'K', 'BIC').opts(size=10, width=500)

El BIC disminuye mientras más simple es el modelo y mientras mayor sea su log-verosimilitud. El mejor modelo es aquel que tiene menor BIC.

En este caso el mínimo BIC se obtiene con \(K=3\). Utilicemos este valor e inspeccionemos los resultados

  • Utilizamos el método predict para obtener las predicciones duras (asignaciones)

  • Utilizamos el método predict_proba para obtener las predicciones suaves (probabilidades)

gmm = GaussianMixture(n_components=3, covariance_type=cov, n_init=3)
gmm.fit(data)
hatr = gmm.predict(data)

The hard assignments:

Hide code cell source
hv.Scatter((data[:, 0], data[:, 1], hatr), vdims=['y', 'z'],).opts(width=500, color='z', cmap='Category10')

Con las predicciones probabilísticas podemos calcular la entropía de Shannon:

\[ H = -\sum p_i \log p_i \]
p = gmm.predict_proba(data)
H = -np.sum(p*np.log(p), axis=1)
Hide code cell source
hv.Scatter((data[:, 0], data[:, 1], H), vdims=['y', 'z'],).opts(width=500, color='z', cmap='Blues', colorbar=True)

Mientras mayor sea la entropía (color más azul) mayor es la incertidumbre respecto a cual cluster pertenece el dato.

20.5.1. Relación con K-means clustering#

El algoritmo K-means (wiki):

../../_images/k_means_algo.png



El algoritmo K-means es un caso particular de GMM con las siguientes restricciones:

  • Clusters con covarianza esférica: \(\Sigma = \sigma^2 \text{I}\), \(\sigma \in \mathbb{R}_{+}\)

  • Clusters con igual varianza: \(\sigma^2 = \sigma_1^2 = \sigma_2^2 = \ldots = \sigma_K^2\)

  • Prior uniforme para la concentración de los clusters: \(\pi_k = \frac{1}{K}\)

  • Se usan siempre asignaciones duras en lugar de probabilísticas

Con lo anterior la regla de actualización de la responsabilidad se reduce a:

\[\begin{split} \begin{split} \hat r_{ik}^{\text{old}} &= \text{arg} \max_k \frac{ \exp \left(-\frac{1}{2\sigma^2} \|x_i - \mu_k\|^2 \right) }{ \sum_{k=1}^K \exp \left(-\frac{1}{2\sigma^2} \|x_i - \mu_k\|^2 \right) } \\ &= \text{arg} \max_k \log \exp \left(-\frac{1}{2\sigma^2} \|x_i - \mu_k\|^2 \right) \\ &= \text{arg} \min_k \frac{1}{2\sigma^2} \|x_i - \mu_k\|^2 \\ &= \text{arg} \min_k \|x_i - \mu_k\|^2 \end{split} \end{split}\]

Y el único parámetro que es necesario actualizar son las medias de los clusters:

\[ \mu_{k}^{\text{new}} = \frac{ \sum_{i=1}^N \hat r_{ik}^{\text{old}} x_i}{\sum_{i=1}^N \hat r_{ik}^{\text{old}}} \]

Y estas se aplican iterativamente al igual que el con el método de EM.

Aquí es un ejemplo de datos dónde k-means es adecuado (fuente):

../../_images/k_means_good.png



Aquí es un ejemplo de datos dónde k-means no es adecuado, pero GMM es adecuado (fuente):

../../_images/kmeans_vs_gmm.png

Ojo: Los clusters en K-means comparten la misma covarianza, pero los círculos se dibujan según los puntos más lejos que pertenecen a un cluster y no reflejan bien la igualdad de la covarianza.



Advertencia

En general GMM clustering y k-means clustering solo funcionan bien en datos con distribuciones más simples. GMM clustering es mejor que k-means clustering en algunos casos. Veamos la imagen abajo para comparar diferentes métodos (algoritmos) de clustering.

../../_images/clustering_comp.png

(Fuente)

20.6. Uso 2: Estimación de densidad y generación de datos#

Podemos utilizar el método sample para generar datos sintéticos con el modelo ajustado antes.

hatdata = gmm.sample(n_samples=10000)
hatdata
(array([[ 2.06340828, -1.07007695],
        [ 2.86170661,  0.10709019],
        [ 3.25894392,  0.08142743],
        ...,
        [-1.75787252,  1.93051828],
        [-2.03886209,  1.51653256],
        [-0.24115727, -0.03886263]]),
 array([0, 0, 0, ..., 2, 2, 2]))
Hide code cell source
hv.Scatter((hatdata[0][:, 0], hatdata[0][:, 1])).opts(width=500)



Ahora, veamos otro ejemplo en el que los datos no son generados por la mezcla de gaussiana. Aquí veamos cómo usamos GMM para la estimación de densidad y generación de datos. (Basado en el fuente)

Hide code cell source
color_iter = itertools.cycle(["navy", "c", "cornflowerblue", "gold", "darkorange", "blue", "red", "black", "orange", "grey"])

def plot_data(X, index, title):
    plt.subplot(1, 1, 1 + index)
    plt.scatter(X[:, 0], X[:, 1], s=10)
    plt.xlim(-6.0, 4.0 * np.pi - 6.0)
    plt.ylim(-5.0, 5.0)
    plt.title(title)
    plt.xlabel("X1")
    plt.ylabel("X2")

def plot_results(X, Y, means, covariances, index, title):
    splot = plt.subplot(1, 1, 1 + index)
    for i, (mean, covar, color) in enumerate(zip(means, covariances, color_iter)):
        v, w = linalg.eigh(covar)
        v = 2.0 * np.sqrt(2.0) * np.sqrt(v)
        u = w[0] / linalg.norm(w[0])
        # as the DP will not use every component it has access to
        # unless it needs it, we shouldn't plot the redundant
        # components.
        if not np.any(Y == i):
            continue
        plt.scatter(X[Y == i, 0], X[Y == i, 1], s=10, color=color)

        # Plot an ellipse to show the Gaussian component
        angle = np.arctan(u[1] / u[0])
        angle = 180.0 * angle / np.pi  # convert to degrees
        ell = mpl.patches.Ellipse(mean, v[0], v[1], angle=180.0 + angle, color=color)
        ell.set_clip_box(splot.bbox)
        ell.set_alpha(0.5)
        splot.add_artist(ell)

    plt.xlim(-6.0, 4.0 * np.pi - 6.0)
    plt.ylim(-5.0, 5.0)
    plt.title(title)
    plt.xlabel("X1")
    plt.ylabel("X2")

def plot_samples(X, Y, n_components, index, title):
    plt.subplot(1, 1, index)
    for i, color in zip(range(n_components), color_iter):
        if not np.any(Y == i):
            continue
        plt.scatter(X[Y == i, 0], X[Y == i, 1], s=10, color=color)

    plt.xlim(-6.0, 4.0 * np.pi - 6.0)
    plt.ylim(-5.0, 5.0)
    plt.title(title)
    plt.xlabel("X1")
    plt.ylabel("X2")
    
# ======= Generar los datos (N=100, D=2) para entrenar GMM ========
# Parameters
n_samples = 100

# Generate random sample following a sine curve
np.random.seed(0)
X = np.zeros((n_samples, 2))
step = 4.0 * np.pi / n_samples

for i in range(X.shape[0]):
    x = i * step - 6.0
    X[i, 0] = x + np.random.normal(0, 0.1)
    X[i, 1] = 3.0 * (np.sin(x) + np.random.normal(0, 0.2))

plt.figure(figsize=(10, 2.5))
plt.subplots_adjust(bottom=0.04, top=0.95, hspace=0.2, wspace=0.05, left=0.03, right=0.97)
plot_data(X, 0, "Datos de entrada (N=100, D=2)")
../../_images/c62770664c23f40206411194e0615b1592e4d448d1211c429616837f7c4d3e36.png
# ========== Entrenar GMM =============
gmm = GaussianMixture(n_components=10, covariance_type="full", n_init=3, random_state=0).fit(X)

plt.figure(figsize=(10, 2.5))
plt.subplots_adjust(bottom=0.04, top=0.95, hspace=0.2, wspace=0.05, left=0.03, right=0.97)
plot_results(X, gmm.predict(X), gmm.means_, gmm.covariances_, 0, "Estimación de densidad por GMM (K=10)")
../../_images/382707c85f51906000954bea8c73051b0a275870275f8e95067caabde67e5bbe.png
# ========== Generar datos por GMM =============
X_s, y_s = gmm.sample(100)

plt.figure(figsize=(10, 2.5))
plt.subplots_adjust(bottom=0.04, top=0.95, hspace=0.2, wspace=0.05, left=0.03, right=0.97)
plot_samples(X_s, y_s, gmm.n_components,1, "Datos generados por GMM (N=100)")
../../_images/392a0277ec410d344fa70414ca3b1ecb9ce331a0c58e15e74f3e5bdd24a9ab0c.png

20.6.1. Relación con Kernel density estimation (KDE)#

Recuerda que la distribución unidimensioanl aproximada con KDE para una muestra \(\{x_i\}_{i=1,\ldots, N}\) es

\[ \hat f_h(x) = \frac{1}{Nh} \sum_{i=1}^N \kappa \left ( \frac{x - x_i}{h} \right) \]

donde \(h\) es el ancho de banda de kernel (kernel bandwidth) y \(\kappa(u)\) es una función de kernel.

Un kernel muy utilizado es kernel Gaussiano:

\[ \kappa(t) = \frac{1}{\sqrt{2\pi}} \exp \left ( - \frac{t^2}{2} \right), \]

Por lo tanto tenemos:

\[\begin{split}\begin{split} \hat f_h(x) &= \frac{1}{Nh} \sum_{i=1}^N \frac{1}{\sqrt{2\pi}} \exp \left( - \frac{(\frac{x - x_i}{h})^2}{2} \right) \\ &= \frac{1}{N} \sum_{i=1}^N \frac{1}{h\sqrt{2\pi}} \exp \left( - \frac{(x - x_i)^2}{2h^2} \right) \end{split}\end{split}\]

Y cómo es la densidad de distribución normal?

\[\begin{split}\begin{split} f(x) = \frac{1}{\sigma\sqrt{2\pi}} \exp(-\frac{(x-\mu)^2}{2\sigma^2})\\ \end{split}\end{split}\]

Asi que tenemos:

\[\begin{split} \hat f_h(x) &= \sum_{i=1}^N \frac{1}{N} \mathcal{N}(x; x_i, h) \end{split}\]

Cómo se relaciona con GMM?

La densidad de un punto de dato se expresa así con GMM (considerando el caso unidimensional):

\[f(x ; \theta) = \sum_{k=1}^K \pi_k \mathcal{N}(x; \mu_k, \sigma_k)\]

Cuál es la ventaja y desvantaja de cada uno?



Material adicional y referencias

Puede profundizar revisando el capítulo 11 del libro “Machine Learning: A Probabilistic Perspective”, de Kevin Murphy, y/o el capítulo 9 “Mixture Models and EM” del libro “Pattern Recognition and Machine Learning” de Christopher Bishop.